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I. Summary 

The recently developed second-order explicit shock-capturing methods of the van Leer, 
Harten, and Yee types, in conjunction with the generalized flux- vector splittings of Vinokur 
and Montagne, and a generalized Roe’s approximate Riemann solver of Vinokur for a real gas 
are studied. The comparisons are made on different one-dimensional Riemann (shock-tube) 
problems for equilibrium air with various ranges of Mach numbers, densities, and pressures. 
Six different Riemann problems are considered. These tests provide a check on the validity of 
the generalized formulas, since theoretical prediction of their properties appears to be difficult 
because of the non-analytic form of the state equation. The numerical results in the super- 
sonic and low-hypersonic regimes indicate that these approaches produce good shock-capturing 
capability and that the shock resolution is only slightly affected by the state equation of equi- 
librium air. The difference in shock resolution between the various methods varies slightly 
from one Riemann problem to another, but the overall accuracy is very similar. For the one- 
dimensional case, the relative efficiency in terms of operation-count for the different methods 
is within 30%. The main difference between the methods lies in their versatility in being 
extended to multidimensional problems with efficient implicit solution procedures. 


II. Introduction 

Several newly developed high-resolution shock-capturing methods [1-5] have been shown 
to be applicable to many multidimensional fluid dynamics problems for a perfect gas. For 
problems containing moderate to fairly strong shocks, these methods produce highly accurate 
solutions near discontinuities [6-10]. These numerical methods belong to the class of total 
variation diminishing (TVD) schemes. 


fAn abbreviated version will appear in the proceedings of the 7th GAMM Conference on Numerical Methods 
in Fluid Mechanics, September 9-11, 1987, Louvain-la-Neuve, Belgium. 
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There exists many ways to achieve higher-order spatial accuracy and at the same time have 
TVD-type properties. Here, two ways are considered. The first way is due to Harten [3], 
Roe [4] and Yee [5]; the second way is due to van Leer [1] (it is sometimes referred to as the 
MUSCL approach). Hereinafter, we will refer to the first way as the non-MUSCL approach. To 
annotate the various algorithms, these nomenclatures will be used. It is emphasized that the 
basic high-resolution shock-capturing methods for hyperbolic conservation laws are developed 
for nonlinear scalar hyperbolic conservation laws. Extension of scalar methods to nonlinear 
systems is accomplished by assuming certain physical models or by local linearization. The 
mathematical foundation relies mainly on the scalar case. There is no identical theory for non- 
linear systems or for the multidimensional counterpart. These schemes are formally extended 
to one- or higher-dimensional nonlinear systems of hyperbolic conservation laws via the so 
called Riemann solvers and are evaluated by numerical experiments. 

There exist three popular ways of extending scalar schemes to nonlinear systems via the 
Riemann-solver approaches: the exact Riemann solvers [11], the approximate Riemann solvers 
[12-15], and the flux-vector splitting techniques [16-18]. However, these Riemann solvers as 
originally developed, are only valid for a perfect gas. An exact Riemann solver [19], approx- 
imate Riemann solver [20-22] and flux-vector splitting approaches [23-25] have recently been 
generalized to a real gas (see reference [7] for details). The objective of this paper is to in- 
vestigate the applicability of second-order explicit shock-capturing methods of the van Leer, 
Harten, Yee types [1,3-5] in conjunction with recently developed generalized flux-vector split- 
tings of Vinokur and Montagne [25] and a generalized Roe’s approximate Riemann solver of 
Vinokur [22] for a real gas. 

The combination of the three Riemann solvers and of the differencing algorithms considered 
above yields five different schemes: a symmetric non-MUSCL scheme [5], an upwind non- 
MUSCL scheme [7], and three MUSCL-type schemes, depending on the Riemann solvers. The 
present study provides a check on the validity of the generalized formulas, since theoretical 
prediction of their properties appears to be difficult because of the non-analytic form of the 
state equation. The values of the state equation are obtained using a curve-fit procedure given 
in reference [26]. Comparisons are made on the accuracy and robustness of the methods. The 
six test cases chosen here are intended to highlight the effect of the high ratios in pressure or 
density related to shocks, and the effect of departure from perfect gas in the state equation. 

A brief description of the Riemann solvers is given in the next section, and the numerical 
algorithms used in conjunction with the various Riemann solvers are then presented. The test 
cases and numerical results follow. 


III. Riemann Solvers 

The conservation laws for the one-dimensional Euler equations can be written in the form 


d JL , jW) n 

dt ^ dx 


where the column vectors U and F(U) take the form 


(3.1a) 
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(3.1b) 


u = 

p 

m 

, F = 

pu 

mu + p 


e 


eu + pu 


Here p is the density, m — pu is the momentum per unit volume, p is the pressure, e = p(e+ 1 it 2 ) 
is the total internal energy per unit volume, and e is the specific internal energy. 

Many Riemann solvers make use of the eigenvalues and eigenvectors of the Jacobian matrix 
A — dF/dU . For a general gas, one therefore requires the thermodynamic derivatives of p . 
In terms of the internal energy per unit volume e = pe, the thermodynamic derivatives can be 
defined as 


f dp\ 

x= ,«p)r 1 


«=(£) 


(32) 


If h = e + p/ p is the specific enthalpy, one can obtain for the speed of sound c the relation 


c 2 = x + K h- 


(3.3) 


For a perfect gas, x = 0, and k = (7 — 1). 
The Jacobian matrix A takes the form 


A = 


0 1 0 
X - (2 - k ) u 2 /2 (2 - /c)u k 

(x + ku 2 /2 - H)u H — ku 2 (1 + k)u 


(3.4) 


where H — h + u 2 /2 is the total enthalpy. The three eigenvalues of A are 

a 1 = u - c, a 2 = u, and a 3 = u + c. 

The corresponding right-eigenvector matrix is 


(3.5) 
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where = ~ 

[ku 2 /2 

+ x] an< l 

b% 

= 

k/c 2 . 

In 

order 

to relate the 

variables p and 

c to the 


independent variables p and c, it is customary to introduce the nondimensional thermodynamic 
variables 


7=1 + “, 

pc 


7 = 


PC 


(3.8) 
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Since there is an arbitrary additive constant in the definition of e, 7 does not have a well 
defined value. In practice, it is introduced only when the arbitrary constant has been chosen 
so that € — ► 0 when p — ► 0. The expressions for flux-vector splittings derived below will be 
given in term of a general e. Only in the numerical results will 7 be employed. For a perfect 
gas the variables 7 and 7 are constant and equal to each other. For an equilibrium real gas 
they are both arbitrary functions of p and e. 


3.1 An Approximate Riemann Solver (Generalized Roe Average [22]) 

Among the various approximate Riemann solvers for a perfect gas [12-15], the most common 
one uses the Roe average [12] because of its simplicity and its ability to satisfy the jump 
conditions. The generalization of the Roe average for a real gas has been carried out by Vinokur 
[22] for the three-dimensional case. A summary of the generalized Roe-average derivation 
for the one-dimensional case is given below. The use of the approximate Riemann solver in 
conjunction with the numerical schemes will be discussed in sections 4.1 and 4.2. 

The flux at a point separating two states U L and U R is based on the eigenvalues and 
eigenvectors of some average A. The optimum choice for A is one satisfying 

A F = AAt/, (3.9) 

where A(*) = (*) fi - (*) L . This choice of A will capture discontinuities exactly. One way of 
obtaining A is to seek an average state t/, such that 

A = A(U). (3.10) 

Such a state is known as a Roe-averaged state. Expressions for a perfect gas were first devised 
by Roe [12]. 

The entries in A depend explicitly on the thermodynamic variables /i, x> and k. Since the 
density is not explicitly required, one would expect the Roe-averaged state to depend on p L and 
p R through their ratio only. It is, therefore, convenient to define the parameter D — p R / p L . 
We first examine the second component of equation (3.9). The average velocity u must be a 
linear combination of u L and u R . Recalling that u L and u R can be independently prescribed, 
we can readily establish the same u as for a perfect gas: 

u L + Du r 

u — 

1 + D 

This definition will satisfy all the terms involving the velocity. Note that u always lies between 
u L and u R . The remaining terms in the equation result in the condition 


(3.11) 


XA p + /cAe = A p. 

This last condition is automatically satisfied for a perfect gas. 


(3.12) 


In order to satisfy the third component of equation (3.9) we also require H to have the same 
form as the perfect-gas version: 

+ (3.13) 

1 + D ' 
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Using the definition of H , equations (3.11) and (3.13) can be combined to define the Roe- 
averaged specific enthalpy as 


h = 


h L + Dh R 

V+D 


D 

2(1 + £>) 2 


(Au) 2 . 


(3.14) 


Note that h could lie outside the range given by h L and h R if Au is sufficiently large. The 
Roe-averaged sound speed is given by equation (3.3) as 


c 2 = x + *h. (3.15) 

For a perfect gas, equations (3.11) through (3.15) are sufficient to define uniquely a 1 , R, and 
R 1 , since \ = 0 and 7c is a given constant. 


For a real gas, equation (3.12) provides only one relation for the variables x and 7c. Since the 
equation of state can be nonconvex, a universally valid unique solution may not be possible. In 
order to gain some insight, we consider the special case in which the states R and L are precisely 
those that satisfy the jump conditions across a discontinuity (either a contact discontinuity or 
a shock wave). Then equations (3.9) through (3.15) are consistent with the exact Riemann 
solver even though x and 7c are not uniquely defined. For a shock wave one obtains 

- h L + D 2 h R 

and 

2 _ Ap 

Ap 

For the special case of a thermally perfect gas, c 2 is a function of h only, and one can readily 
show that the values of h and c 2 given by equations (3.16) can only satisfy this law if the 
function is linear. But this is precisely the definition of a perfect gas. 


(3.16b) 


(3.16a) 


The above analysis makes it clear that for a real gas the values of \ and 7c must be defined 
in terms of the thermodynamic states R and L, and not in terms of the state h. Because of 
the nonconvex nature of a real-gas equation of state, the values of x and n at states R and 
L, or some average of the two, will not satisfy equation (3.12), in general. One way to obtain 
unique values of x and 7c is to project the average state given by R and L onto the straight 
line defined by equation (3.12). Since the value of Ae depends on the arbitrary constant in the 
definition of e, the resulting value of 7c will depend on the choice of this constant. To overcome 
this arbitrariness, one divides equation (3.12) by /c. The straight line for the variables \/k 
and \/ K has a slope given by A p and A p, both of which are uniquely defined by states R 
and L . Since A p and Ap are not dimensionally consistent, one must further introduce a scale 
factor proportional to their ratio in order to have the final answer independent of the choice 
of dimensional units. A natural candidate is the square of an average sound speed. In terms 
of arithmetic averages, and the scale factor c 2 = [( c L ) 2 + {c R ) 2 \/ 2, one obtains the expressions 


1 

7c 


{ \ [(fr + Ap + / ( a p 2 + a p 2 ) 


(3.17a) 


5 


and 


I = { $ i<S + S )Ap + + / (Ap! + ^ Ap!) <3i7b) 

Other expressions can be obtained by using different scalings and taking different types of 
averages. The optimum choice that would be valid for a wide range of conditions will require 
further research and numerical experiments. 

3.2 Generalized Steger and Warming Flux- Vector Splitting [25] 

For subsonic flow, the eigenvalues of A are of mixed sign. In flux-vector splitting methods, 
the flux F is divided into several parts, each of which has a Jacobian matrix whose eigenvalues 
are all of one sign. The flux-vector splitting approach of Steger and Warming [16] made use 
of the assumption that the Euler equations are homogeneous of degree one. This homogeneity 
property is satisfied for a thermally perfect gas. While they only gave the results for a perfect 
gas, the results for a thermally perfect gas are found in reference [25]. The generalization to an 
arbitrary equilibrium real gas presented here makes direct use of the eigenvalues of A , and is 
an extension of the work of Sanders and Prendergast [18] and Montagne [22]. The basic idea is 
the observation that the eigenvalues are actually three velocities. We can associate with each 
eigenvalue a*, l = 1,2,3, a fictitious stream with velocity a*, and some unknown density p l and 
specific internal energy e l . Each stream can then be characterized by a column vector 

(3.18) 

where m l ~ p l a l and e l = p l [e l + |(a z ) 2 ]. The fictitious flux associated with each stream is 
assumed to be convective only, namely, 

F l = a l U l . (3.19) 

The six unknowns are determined from the conditions 

3 3 

U = J2 ul and F = ^2f 1 . (3.20) 

i= 1 1=1 

Since the second component of U and the first component of F give the same equation, we 
are left with one degree of freedom. From the first two components of U and F, one readily 
obtains 

p 1 = P 3 = ^ and p 2 - p( 1 - i). (3.21a) 

The third components of U and F result in the relations 

f 1= e 3 and . (3.21b) 


W = 


m 
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The relations in equations (3.21a) and (3.21b) have the same form as those for the Steger- 
Warming solution for a thermally perfect gas. In choosing the manner in which the remaining 
degree of freedom is expressed, we are guided by the fact that e is not an absolute quantity. 
Consequently, we express the € % as 


(1 - a)c 2 

7(7 - 1) ’ 


and 


e = e 


c + 


1 - a 



(3.21c) 


Here the non-dimensional thermodynamic parameter a has been chosen so that a = 0 corre- 
sponds to the Steger- Warming solution for a thermally perfect gas. The final expressions for 
the F l can be written in the form 


and 


F 2 _ pj 7-1) 
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+ 




U~«)« a 

Tf(Tf-l) 




Isil! 
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a' 

(a ') 2 



(3.22a) 


(3.22b) 


for / = 1 or 3. We thus have a one-parameter family of flux-vector splittings, where a can be 
an arbitrary parameter. The total flux can be split according to the sign of the eigenvalues. 
For —c < u < 0 we, therefore, have 

F + = F 5 and F~ = F l + F 2 (3.23a) 

and for 0 < u < c 

F + = F 2 + F 3 and F~ = F 1 . (3.23b) 


For a perfect gas one can show that when — c < u < 0, the determinant of the split-flux 
Jacobian A + = dF + /dU is 


det(A + ) - 


(c + «) 3 

I67 2 



7 


7 J 


(3.24a) 


The determinant is the product of the three eigenvalues, and a necessary condition for to 
have eigenvalues that are all positive is that det(A + ) > 0. It follows from equation (3.24a) 
that we must take a < 7(3 - 7)/ [2(7 - 1)]. For the region 0 < u < c, one can show that the 
minimum of det(A + ) occurs when u approaches zero. This minimum value is 


det(A + ) 



5 — 37 — 2a 


(3 - 7)(7 - 



(3.24b) 
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The condition det(A + ) > 0 requires that a < -y(5 - 37)/ [2(3 - 7)(7 - 1)]. For 1 < 7 < 3, this 
second condition on a automatically satisfies the previous inequality derived for -c < u < 0. 
In order to determine an optimum value of a, we note that in general the three eigenvalues 
of A+ are discontinuous at -c, 0, and c. For a = 0, two of the eigenvalues are continuous 
at u — c, and one can further show that all the eigenvalues are real and non-negative in 
the subsonic region when 7 < 5/3. For a thermally perfect gas with vibrational excitation, 
one can show that all the eigenvalues remain real and non-negative when a = 0. Because of 
the general nature of the variation of 7 and e, it is simplest to set a = 0 for an arbitrary 
equilibrium real gas. The numerical results presented below were actually based on an earlier 
parameterization which was chosen to reduce to the perfect-gas solution. It corresponds to 
the present a = (7 — 7)/ (7 - 1). The cases were also rerun with a = 0, and no noticeable 
difference was observed. 


3.3 Generalized van Leer Flux- Vector Splitting [25] 

In a different approach, van Leer [17] constructed a flux- vector splitting for a perfect gas 
in terms of low-order polynomials of u, which gives continuous eigenvalues at u = 0 and 
u — ±c. The splitting also has the desirable property that one of the eigenvalues of the split 
-flux Jacobians is identically zero. This results in a sharper capture of transonic shocks. An 
extension of van Leer’s splitting for a real gas was derived by Montagne [24], but it is not 
internally consistent. The present formulas represent the most general, consistent extension 
of van Leer’s splitting for a real gas. Because of the arbitrary nature of a real-gas law, the 
condition of one eigenvalue being identically zero cannot be satisfied exactly. 

For \u\ < c , the continuity requirements necessitate a factor (u ± c) 2 in the formulas for 
F ± . The expressions for the first two components of F ± that are given by the lowest-order 
polynomials in u are readily found to be 

fS... = ±£(«±')’ (3.25a) 


and 

FLm = ± c) 2 [2 c ± (7 - l)u] (3.25b) 

For a perfect gas they are identical to those derived by van Leer. The expression for the third 
component, which satisfies continuity and symmetry conditions, can be written most generally 
in the form 


F* =± — 

ene 2(7 2 -1)4 c 


(u ± c) 2 [2c ± (7 — l)til 2 ± y- (u ± c) 
L 4c 




7 ( 7 - 1 ). 


±7 £ (u±c) 2 ( u Tcf, 


(3.25c) 


where the arbitrary parameter 0 has been chosen so that 0 ~ 0 corresponds to van Leer’s 
solution for a perfect gas. We again have a one-parameter family of flux-vector splittings. 

For a perfect gas, the second term in Ff^ t vanishes. Van Leer’s condition of a vanishing 
eigenvalue for A ± requires 0 to vanish also, so that F± nt reduces to one term. One can readily 
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show that the remaining eigenvalues are both of the proper sign for 1 < 7 < 3. For a real gas, 
one requires at least two terms for F e ^ e . (The one-term solution of Montagne [24] actually 
violates the condition F = F+ + F~). Because of the general nature of the variations of 7 
and e, it is impossible to obtain the vanishing eigenvalue condition identically throughout the 
velocity range, for any choice of /?. 

One can demonstrate readily that for a thermally perfect gas, the two-term solution has 
one eigenvalue that is of the wrong sign for the whole subsonic velocity range* Fortunately, 
the magnitude of the offending eigenvalue is extremely small, so that its effect on a numerical 
scheme is not noticeable. In view of this fact, and because of the general nature of the variation 
of 7 and e, it is simplest to set /? = 0 for a real gets. The results presented in section V are 
obtained with this value of /? and seem to show that this approximation is valid for a wide 
range of flow conditions. 


IV. Description of the Numerical Algorithms 

Let U* be the numerical solution at x = j Ax, t = nAt, with Ax the spatial mesh size 
and At the time-step. The second-order in time and second- or third-order in space, explicit 
difference schemes considered here for both the MUSCL and the non-MUSCL approaches can 
be written in the following form 

= ,), (4.1) 

where A = At/ Ax. The vector functions Hj± i are sometimes referred to as numerical flux 
functions. 


4.1 Non-MUSCL Approach 

The numerical flux function Hj± i for a non-MUSCL-type approach for both the upwind 
and symmetric TVD schemes [6] using an approximate Riemann solver can be expressed as 


*** = 5 


F(U,) + F(U J+1 ) + R j+l * j+ , . 


(4.2) 


Here Rj+± is the eigenvector matrix of the Jacobian A evaluated at some symmetric average 
between Uj and 


Second-order Symmetric TVD Scheme : The elements of the + i denoted by + for a 
general second-order symmetric TVD scheme [5] are 


(< + i) S = -A(a' + i )*Q‘ + i - + [a' + i -Q' + .] • 


(4.3a) 


The value a* + i is the characteristic speed a evaluated at some symmetric average between Uj 
and Uj +1 . The function tf> is 



[z' + SD/ 26 , 


\z\ > 6i 

|*i < 6i ' 


(4.3b) 
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Here ip(z) in equation (4.3b) is an entropy correction to \z\ where S\ is a small positive pa- 
rameter (see references [7,27] for a formula for 6x). Since all of the test problems contain 
unsteady shocks, 61 is set to zero in all of the computations. For steady-state problems con- 
taining strong shock waves, a proper control of the size of 61 is very important, especially for 
hypersonic blunt-body flows. See reference [7] for detailed discussion and numerical examples. 
Examples of the limiter function Q l j+x can be expressed as 


Q‘ J+ x = minmod (a' _ i , <*' + 1 ) + minmod (aj + 1 , + § ) - a' + i 
Q' + i =minmod(a'_i,a' + i,a' + | ) 

Q‘ + i =minmod 2 a' _ ±,20^,2^, + a' + § ) . 


(4.3c) 

(4.3d) 

(4.3e) 


The minmod function of a list of arguments is equal to the smallest number in absolute value 
if the list of arguments is of the same sign, or is equal to zero if any arguments are of opposite 
sign. Here + x are elements of 


* = r J+l(Uj+i - Uj). 


(4.4) 


Second-Order Upwind TVD Scheme : The elements of the denoted by (<^ + 1 ) L for a 

second-order upwind TVD scheme, originally developed by Harten, and later modified and 
generalized by Yee [7,28], are 


^ 0 


(4.5a) 


(4.5b) 


(4> l J+ xf = )(»}+! + 9j) - ^(a' + i+l' + iK+i- 

The function a(z) = |[t/>(z) — Az 2 ] and 

, , , . i (»;■+!- *})/°i+j 

-».■+!= '(Vi) 1 0 

Examples of the limiter function g l 3 can be expressed as 
9j = minmod (a' _ x,<x l }+ x) 

gl 3 = ( a y+* a i-i + K+* a i-*l) / ( a i+i +a !-i) 

g J = { a i-i + 8 *\ + [K-i) 2 + ^]}y 

g l j - S •max|o,min(2|a' + i|,5 • a'_ i ),min(|a' + ± |, 2S • a'_±) j; S = sgn(a‘ + i).(4.5f) 
In limiter (4.5e) 5 2 is a small parameter (see reference [1]). 


K + *) , + K -*)*+2 * 


(4.5c) 

(4.5d) 

(4.5e) 
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4.2 MUSCL Approach 

A/ 

MUSCL Approach Using an Approximate Riemann Solver : The numerical flux function 1 
for a MUSCL-type approach of an upwind scheme as described in Yee [29] using an approximate 
Riemann solver can be expressed as 





L 

3 + 2 




(4.6a) 


where the elements of 



a are 
2 


l I -M I 

- — i a‘ i 
J • 2 J'* 2 


4+* 


S J+ * =«7 + , i (i f f + * - 


i+ ; 


J + 5 


(4.6b) 

(4.6c) 


Here a 1 . A are the eigenvalues and f?, . a. is the eigenvector matrix of A, evaluated using an 

3 * 2 ' 2 

approximate Riemann solver between the two states L and that is, 


S; + * (4.6d) 

j)- <<•«•) 


The variables l/* A and Uf A are obtained by interpolating neighboring state variables, 
using limiters. However, there are options in choosing the types of dependent variables in 
applying limiters for system cases. Namely, for the Euler equations, one can impose limiters 
on the conservative, primitive, or characteristic variables. For the numerical tests only one 
option has been used, imposing the limiters on the set of variables W = (p,u ,p(). Then, if 
we let P and F -1 denote transformation operators such that U — PW and W = P~ l U, the 
vectors U R . and U L A for a second-order in time, second-order in space MUSCL approach 
can be defined as 




i+ 1 


n+i v 

+ 






(4.7a) 

(4.7b) 


Here gj is defined as in equations (4.5c)-(4.5f), except the arguments will be (W" +1 - W") and 
[Wp - Wp_ x ), with W ” + * = p- 1 U” + = , where 


t/ n+ " = J7»- - 
J j 2 


r (/’(»? + jw) - f (m? - 5?i)) 


(4.7c) 
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/N/ 

MUSCL Approach Using Flux-Vector Splittings'. The numerical flux + i for either 

flux-vector splitting can be expressed as 

( 4 . 8 ) 

where F ± (U L ' I l) are evaluated using either the generalized Steger- Warming flux-vector split- 
ting or the generalized van Leer flux- vector splitting. The vectors U R , and U L . are the 

Jr 2 ' o 

same as in equations (4.7). 


V. Test Cases 


Six test cases will be presented for various combinations of Mach numbers, densities and 
pressures. They are all one-dimensional Riemann problems. The left- and right-hand-side 
states of the initial conditions are tabulated in table 1. The cases have been ordered in the 
direction of increasing maximum Mach numbers encountered in the flow. 


Table 1. Initial conditions for the test cases. 


Case 

State 

Density, 
kg/m 3 

Pressure, 

N/m 2 

Temp., 

K 

Energy, 

(m/sec) 2 

Velocity, 

m/sec 

Mach No. 

Case A 

Left 

0.0660 

9.84 10 4 

4390 

7.22 10 6 

0 

0.0 


Right 

0.0300 

1.50 10 4 

1378 

1.44 10 6 

0 

0.0 

Case B 

Left 

1.4000 

9.88 10 5 

2438 

2.22 10 6 

0 

0.0 


Right 

0.1400 

9.93 10 3 

2452 

2.24 10 6 

0 

0.0 

Case C 

Left 

1.2900 

1.00 10 5 

272 

1.95 10 s 

0 

0.0 


Right 

0.0129 

1.00 10 4 

2627 

2.75 10 6 

0 

0.0 

Case D 

Left 

1.0000 

6.50 10 5 

2242 

2.00 10 6 

0 

0.0 


Right 

0.0100 

1.00 10 s 

346 

2.50 10 5 

0 

0.0 

Case E 

Left 

0.0100 

5.73 10 2 

199 

1.44 10 5 

2200 

7.8 


Right 

0.1400 

2.23 10 4 

546 

4.00 10 5 

0 

0.0 

Case F 

Left 

0.0100 

5.73 10 2 

199 

1.44 10 5 

4100 

14.6 


Right 

0.0100 

5.73 10 2 

199 

1.44 10 5 

-4000 

-14.5 


The thermodynamic properties of equilibrium air are obtained from the curve fits of Srini- 
vasan et al. [26]. These curve fits give analytic expressions for 7 in several ranges of density 
and internal energy. The values of 7 are then calculated from the derivatives of these analyt- 
ical expressions. The numerical solution is compared with an “exact solution” computed by 
solving the Rankine-Hugoniot jump conditions and integrating numerically the characteristic 
equations in the expansion fan. Figure 1 shows the distribution of Mach number and of the two 
quantities 7 and 7 defined in equation (3.8). Not only the changes in their values but also the 
differences between them are an indication of departure from the perfect gas case. These dif- 
ferences do not necessarily occur at very high temperatures, but at intermediate temperatures 
when the vibration is excited and when the dissociation reactions start. 

The six cases have also been tested for a perfect gas (7 = 1.4), with the same initial conditions 
of density and internal energy. Indeed, though they correspond to significant variations of the 
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the thermodynamic properties of the gets, the conditions of test cases A, B, and D are not very 
severe. On the other hand, test case C corresponds to a large jump of internal energy through 
the contact discontinuity, and cases E and F present strong shocks at different Mach numbers. 
Since real-gas effects occur in general in flows with large discontinuities, it is important to 
identify this last effect in order to estimate the performances of the methods for the real gases. 


VI. Numerical Results 

Five different second-order explicit schemes have been tested for the cases described in table 
1, namely, 

a. The symmetric TVD scheme (4.3), non-MUSCL 

b. The upwind TVD scheme (4.5), non-MUSCL 

c. The upwind TVD scheme (4.6), MUSCL 

d. The generalized van Leer flux-vector splitting (4.8), MUSCL 

e. The generalized Steger- Warming flux- vector splitting (4.8), MUSCL 

Schemes (a) and (b) follow the non-MUSCL approach, schemes (c)-(e) follow the MUSCL 
approach. For the test cases presented here, the same approximate Riemann solver is used in 
schemes (a)-(c). The distribution of six variables 7, 7, u/c, e, u, and p are shown together 
because the differences between the schemes do not always occur for the same variable. After 
a presentation of the results obtained with these five schemes, some issues concerning the 
comparison of the differencing methods and of the Riemann solvers are addressed. 

6,1 Numerical Experiments 

The numerical experiments were done with a uniform grid. The number of discretization 
points is 141 in cases A-D (table 1), and only 81 for cases E and F because in the latter cases the 
expansion fan is replaced by a shock. The origin x = 0 corresponds to the initial discontinuity 
between the left and right states. In the grid, this origin corresponds to the middle of the 
interval number 70 for cases A-D, 16 for case E, and 35 for case F. The time for stopping the 
computation has been chosen for each case in order to use the full computational domain. In 
the case of equilibrium air, for a Ax normalized to 0.1 m, these stopping times, are t = 3.5 
msec for case A, t = 3.0 msec for case B, t = 5.0 msec for case C, t = 3.2 msec for case D, 
t ~ 7.0 msec for case E and t = 6.0 msec for case F. These times have not been changed for 
the perfect gas except that t = 4.0 msec for case F because the differences in 7 dramatically 
change the shock motion. 

The whole set of results presented has been obtained by using the same limiter function for 
each scheme. Scheme (a) uses the limiter (4.3e). Schemes (b)-(e) use the limiter (4.5e). For 
the numerical values used in all the test cases presented, a value of 10 -7 for an appropriately 
normalized parameter £2 was found to be adequate. Results using scheme (b) with the limiter 
(4.5d) on the nonlinear fields and limiter (4.5e) on the linear field can be found in references 
[5,7], 

The results presented for the approximate Riemann solver have been obtained with two 
different averages. The results for cases A, D, and F used the generalized Roe average (section 
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3.1), with one modification. The projection was done on the straight line defined by equation 
(3.12), and the scale factor was chosen as h. Thus equations (3.17a) and (3.17b) were replaced 

by 


= = | ^ [— — -= - ^ — At - (k l + K R )hAp] A( + Ap/iAp j j {At 2 + ( hAp ) 2 ) (6.1a) 


and 


k — |ApA?- ^ Ac - {k l + K R )hAp\hAp^ j (Ae 8 + (hAp) 2 ). (6.1b) 

The results for the cases B, C, and E used modifications of the average of Carofano [21]: 


2 - 0X& 


5 + (ffll) 

h h L + h* ' 

k = 0.5 (k l + k r ) 

X = c 2 - nh. 


(6.2a) 

(6.2b) 

(6.2c) 


This average recovers Roe’s average for a perfect gas but does not satisfy Roe’s property ”U” 
for a real gas. 

The time-step limit is expressed in term of a CFL number that has been defined as the 
maximum value of (|u| + c) At/Ax. For those cases in which the maximum value occurs in a 
subsonic region, such as in test case C, it was found that the flux-vector splitting requires a 
more stringent CFL limit. When |u/c| is less than one, the value of (|u| + c)At/Ax for the 
generalized van Leer flux-vector splitting had to be multiplied by the factor [17] 

ki = (|u| + c)( 7 + Z) j [2-y + |u/c|(3 - 7)] (6.3) 

whereas for the Steger- Warming flux-vector splitting a less stringent modification was found 
adequate, using a factor 

k 2 = lj{0.8 + 0.2\u/c\). (6.4) 

This CFL number has been set equal to 0.9 except for case C. In case C it has been fixed at 
0.5 for the upwind non-MUSCL and MUSCL schemes, and at 0.9 for the symmetric scheme. 
We briefly describe here the results obtained for each case. 

Test case A : The variations of 7 and 7 are important. The plots for test case A (figs. 2, 3) 
show that for the intermediate states of the numerical solution at the contact discontinuity, the 
variations of 7 are not even monotonic. Actually, in the Mollier diagram, 7 is not a monotonic 
function of energy, and in that particular case, the two end states of the contact discontinuity 
are located on either side a minimum for 7. Even though the distribution of such quantities 
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as density or energy is monotonic, an oscillation necessarily occurs for 7. This behavior may 
be the cause of a small wiggle for the velocity distribution at the contact discontinuity, which 
does not appear for the perfect gas. Otherwise, the comparison between the perfect gas and 
the real gas indicates no deficiency in the generalization techniques. The differences between 
the schemes are very small for this problem. They all capture the shock within three points. 
The best results are obtained with the schemes (b), (c) and (d). No significant difference 
appears between the symmetric scheme (a) and the upwind schemes, either non-MUSCL (b) 
or MUSCL (c). 

Test case B : The intensity of the shock in case B is greater than in case A (figs. 4, 5). The 
capture of the shock is slightly sharper with the schemes (a)-(c) using the approximate Riemann 
solver. In this case, no particular problem can be detected at the contact discontinuity. The 
symmetric scheme (a), however, captures it within more points than the upwind schemes. 
Looking in particular at the velocity distribution, the most regular solution is given by the 
scheme (c) in this case. Again, the behavior of all the numerical solutions are quite similar 
using either the perfect gas or the real gas. 

Test case C : Case C is a difficult test for all schemes, either for the real (fig. 6) gas or perfect 
gas (fig. 7). There are noticeable errors on the levels of the velocities and on the position 
of the shock. The errors are a little greater for the real gas. All the schemes have problems 
at the contact discontinuity for the real gas. Especially for scheme (e), using the generalized 
Steger- Warming flux splitting, the value 0.9 of the CFL number without the correction yields 
wiggles that propagate everywhere on the solution. For the perfect gas, only schemes (d) and 
(e) using the flux-vector splitting technique have a problem. This may be related to the fact 
that unlike for the five others cases for the perfect gas, the absolute value of the Mach number 
crosses 1.0 at the contact discontinuity. At the shock also, the schemes have some problems. 
The resolution is very sensitive to the time-step, especially for scheme (b). This was the reason 
for choosing the CFL value of 0.5 for the upwind schemes and 0.9 for the symmetric scheme, 
which seems to give the best results at the shock. As the same problem occurs both for the 
real gas and the perfect gas, this case shows the limits of the numerical methods even more 
than those of their extensions to a real gas. 

Test case D : Case D (fig. 8, 9) is characterized by a large transonic rarefaction wave. The 

discontinuity, which occurs in the middle for most of the variables on the real gas case, is 
also present on the distribution of 7 and 7 for the exact solution. The tables from reference 
[26] fit the state equation with continuous functions on a set of regions in density and energy 
which covers the physical range of interest. However, continuity is only loosely ensured at 
the boundaries of these regions. The discontinuity observed here corresponds to a density 
boundary. Here, the discontinuities in the values given by the curve fits for the pressure and 
its derivatives is approximately 5%. In that case, only scheme (e) using the generalized Steger- 
Warming flux-vector splitting recovers a smooth solution at that point. As in case C, but 
with less intensity, an error on the levels of energy and velocity appears, either for the real 
gas (fig. 8) or the perfect gas (fig. 9). This error is less important with the non-MUSCL 
schemes (a) and (b). The results are shown for a value of the CFL number equal to 0.9, which 
corresponds to the sharpest resolution of the shock. At a lower value of the CFL number (0.5), 
the energy overshoot at the contact discontinuity is almost completly released for the non- 
MUSCL schemes (a) and (b). On the contrary, this does not help with the MUSCL schemes 
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(c)-(e). 


Test case E : In case E (fig. 10, 11), the strong shock on the left and the contact discontinuity 
correspond to large jumps on the thermodynamic variables. The shock is almost stationary 
and is captured within the same number of points with all the schemes. However, the lowest 
oscillation just behind the shock is achieved by scheme (d) using the generalized van Leer 
flux-vector splitting. The contact discontinuity is best captured with the symmetric scheme 
(a). All the results present some wiggles on the energy and the velocity in the region between 
the shocks. However, the levels are straighter with schemes (a) and (e). For schemes (c)-(e) 
the choice of the limiter function is more critical. They blow up when using limiter (4.5d) or 
(4.5f). Again, the behavior of the schemes is very similar for the real gas and the perfect gas. 
However, for the perfect gas, the wiggles are slightly smaller. 

Test case F : The situation for case F (figs. 12, 13) is almost the same as for case E. The 

approximate Riemann solver captures the shock within one point less than the flux-vector 
splittings. The extra point in the shock given by flux-vector splittings (d) and (e), yields a 
value of 7 that does not lie between the limits on each side. As in case A, this is due to the 
nature of the state equation. The bump observed in the distribution of energy and density is at 
the location of the contact discontinuity. It is mainly a result of errors introduced at the very 
first time-steps, when the whole phenomenon is represented, without enough accuracy, on a 
few space intervals (tests on the same case at lower Mach number show that the error reduces 
with the value of the Mach number). Some schemes seem to damp this error more effectively, 
for example, scheme (e) using the Steger- Warming flux-vector splitting and symmetric TVD 
scheme (a). Moreover, the solutions obtained with symmetric TVD scheme (a) and upwind 
scheme (b) non-MUSCL schemes are very similar since this particular example contains shock 
waves only. Because of the slight asymmetry introduced in the velocities of the initial states, 
the velocity of the intermediate state of the Riemann problem is not exactly zero. However, the 
strenghts of the two shocks are the same and their numerical solution should be symmetrical. 
As can be seen from the density plots, this is not exactly true for schemes (d) and (e) using 
the flux-vector splittings. Moreover, using schemes (d) and (e), a steady solution is formed 
within a few time-steps, with the two initial states separated at one point by a state that has 
a nearly zero velocity. This seems to be due to the fact that the states U R x and U L _ t defined 

by equations (4.7) around this point remain supersonic. In order to overcome this difficulty, 
it has been necessary to impose a zero value on the slope of the velocity at the intermediate 
point during the first time-steps. Otherwise, the effects of limiters and of the state equation 
are similar to case E. 

6.2 Comparison of the Differencing Techniques 

Parts (a)-(c) of figures 2-13 provide comparisons of the symmetric TVD scheme, the upwind 
TVD scheme, and the MUSCL scheme with the same approximate Riemann solver for cases 
A through F (table 1 ). On the basis of the numerical tests, some general observations can be 
made about the numerical resolution and the computational efficiency of the methods. Since 
it was not possible to display the results obtained with all the limiter functions ((4.3c)-(4.3e); 
(4.5c)-(4.5f)) introduced for the different schemes, some additional comments are given about 
their use on the test cases. 
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Numerical Resolution : The results obtained with the differencing techniques rely on the va- 
lidity of the Riemann solver that they use. From this point of view the real-gas case is not 
theoretically different from the perfect gas. Then the conclusions found here for the accuracy 
are mostly related to the flow conditions in Mach number and intensity of the discontinuities. 
For relatively gentle cases like A and B, the techniques are hardly distinguishable. Greater 
differences appear in the representation of the shock between the non-MUSCL and the MUSCL 
approach in case C and also in cases E and F. Particularly in case C, the non-MUSCL approach 
appears to be slightly more sensitive to the time-step. For the treatment of the contact discon- 
tinuities the non-MUSCL approach seems to be slightly better, especially in case D. But the 
major differences are between the symmetric scheme and the two upwind schemes. Although 
the upwind schemes give better solutions for the rapidly unsteady cases like A-D, especially 
for contact discontinuities, the situation is reversed in cases E and F where the motion of the 
shock is slow. The levels are more flat with the symmetric scheme; moreover, the capture of 
the contact discontinuity is sharper with the symmetric scheme in case E. Finally, cases C and 
D seem to be at the accuracy limit of the methods, within the number of points used here for 
the discretization of the problems. 

Influence of the Limiter Functions: The limiter function is an important parameter for the 
TVD schemes. The choices used here correspond to what we considered as a good compromise 
between stability and accuracy on the whole set of test cases. In particular for the cases A, B 
and D, the treatment of contact discontinuities and expansion waves can be improved by using 
limiters (4.5d) or (4.5f) instead of (4.5e) for the upwind schemes (b)-(e). However, when using 
limiter (4.5d) or (4.5f) in cases E and F, some instabilities have been found with the MUSCL 
schemes near the shocks for high Mach numbers. Similarly, limiter (4.5f) with either upwind 
scheme yields an overshoot at the contact discontinuity for case C. The different tests seem to 
show that the influence of the limiter is the same as for a perfect gas. Detailed discussions of 
this subject can be found in [5,7] for the non-MUSCL schemes, and in [1,9] for the MUSCL 
scheme. 

Computational Effort: the use of a curve fit for the computation of thermodynamic variables 
like the pressure and sound speed from the conservative variables can be considerably more 
time consuming than the simple analytic formulas existing for a perfect gas. The non-MUSCL 
approach evaluates these properties once at the center of each cell. On the other hand, the 
need for computing the formulas of the Riemann solvers at each interface requires two such 
evaluations per cell and per space dimension for the MUSCL approach. For the one-dimensional 
tests, this difference was not very critical. For multidimensional calculation, the differences 
can be larger and depend greatly on the efficiency of the software package used for the curve 
fits. 


6.3 Comparison of the Riemann Solvers 

Parts (c)-(d) of figures 2-13 compare the three Riemann solvers for the MUSCL scheme. In 
this part, we are mostly interested in checking the validity of the Riemann solvers. For this 
purpose, we analyze here the effects of the two different phenomena occurring in these flows, 
the large discontinuities and the real gas state equation. 

Influence of the Mach Number and of Large Discontinuities: Since these methods have been 
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used for a long time with a perfect gas, remarks will be made only about the treatment of large 
discontinuities which are found especially in cases C, E and F. All three methods seem to be 
able to handle shocks corresponding to fairly high Mach numbers. However, the wiggles that 
appear in the subsonic region in cases E and F tend to increase with the Mach number, for the 
perfect and real gases. Form this point of view, the flux-vector splitting methods which are 
more dissipative show slightly fewer of these wiggles. On the other hand, case C seems to show 
some limits for these Riemann solvers in transient flows with large discontinuities. In general, 
the approximate Riemann solver is more accurate for the contact discontinuities. Especially 
in case C, it gives a better resolution where flux-vector splitting methods have trouble at the 
contact discontinuity, even with the perfect gas. 

Real Gas Effects : In the test conditions used, the comparisons of the different Riemann 

solvers is not changed very much by the use of equilibrium air. The advantages of the gener- 
alized van Leer flux-vector splitting over the generalized Steger- Warming flux-vector splitting 
remain the same as for the perfect gas in general. The three generalization techniques involve 
a free parameter. Though the choices used in the numerical tests correspond to natural ex- 
tensions of the perfect-gas techniques, a wide variety of possibilities has not been explored. 
The way of determining the best choice for the most accurate and stable results is not clear. 
In the case of the approximate Riemann solver, even the necessity of using an average satis- 
fying Roe’s property ”U” is an open question. The average (6.1) did not yield a noticeable 
improvement on the results obtained with the other one (6.2). Another important fact is that 
flux-vector splittings make use of the sound speed only, whereas approximate Riemann solvers 
of the Roe- type make use of the thermodynamic derivatives \ and * of equation (4.2). These 
thermodynamic derivatives put more stringent requirements on the curve fit that represents 
the thermodynamic properties of the gas. In this regard, the curve fits of Srinivasan et al. 
may be deficient for the approximate Riemann solver as can be seen from figure 5, case D. 
One probably needs more improved curve fits than those of reference [26] before a definite 
conclusion can be drawn on the accuracy of the different Riemann solvers and schemes. From 
this point of view, the generalized Steger- Warming flux-vector splitting, which is more diffusive 
than the two others Riemann solvers, is also less sensitive to a lack of regularity in the curve 
fits. 


VII, Concluding Remarks 

For the purpose of calculations in gas dynamics with an equilibrium real gas, the numerical 
tests show that the simple extensions to a real gas for the flux-vector splitting or the approx- 
imate Riemann solver presented in this paper are valid. The main effect of using a real gas 
equation of state is to exacerbate the problems of the methods for large discontinuities. Test 
case C is an example of such a situation. Similarly, it seems difficult to give a scaling of the 
methods. Depending on the case, each one presents some drawbacks or some advantages. 

Looking at the Riemann solvers in general, the advantages of the generalized van Leer flux- 
vector splitting over the generalized Steger and Warming formulation remain for a real gas, 
with slightly less dissipation at the discontinuities. The approximate Riemann solver gives 
results very similar to the generalized van Leer flux-vector splitting formulation. Looking 
at the numerical schemes, the main difference seems to occur between the upwind and the 
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symmetric TVD schemes. The upwind schemes give better results in strongly transient flows. 
On the other hand the symmetric scheme has a better stability in regions of low Mach number 
and produces shock resolution similar to that of the Steger- Warming flux-vector splitting in 
the MUSCL approach. 

The difference in computational effort required by the three Riemann solvers is small. The 
main difference lies in the MUSCL and non-MUSCL approaches. For one-dimensional prob- 
lems, the operations-count between the non-MUSCL and MUSCL approaches is within 30% 
for a perfect gas. However, due to extra evaluation in the curve fitting between the left and 
right states in an equilibrium real gas for the MUSCL formulation, additional computation is 
required for the MUSCL approach. The effective amount of additional computation increases 
nonlinearly as the spatial dimension increases. None of the differences observed for the ex- 
plicit versions seems to be decisive for the one-dimensional tests, but factors such as stability 
and computational efficiency need further investigation in multidimensional tests, which are in 
progress. Preliminary study shows certain advantages of the approximate Riemann solver over 
the flux-vector splitting approaches. The present results also indicate that the state equation 
does not have a very large effect on the general behavior of these methods for a wide range of 
flow conditions. The ease of extending these methods to multidimensional problems and/or to 
efficient implicit algorithms will be a more important issue for application to hypersonic flow 
calculations. 
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Comparison of the schemes for test case A, for equilibrium air. 

(a) Symmetric TVD, non-MUSCL; (b) Upwind TVD, non-MUSCL; (c) Upwind TVD, MUSCL; 
(d) Generalized van Leer splitting, MUSCL: (e) Generalized Steger and Warming splitting, MUSCL. 

Note: (a)-(c) use the identical approximate Riemann solver. 
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Comparison of the schemes for test case B, for equilibrium air. 

(a) Symmetric TVD, non-MUSCL; (b) Upwind TVD, non-MUSCL; (c) Upwind TVD, MUSCL; 
(d) Generalized van Leer splitting, MUSCL; (e) Generalized Steger and Warming splitting, MUSCL. 
Note: (a)-(c) use the identical approximate Riemann solver. 




Comparison of the schemes for test case B, for the perfect gas. 

(a) Symmetric TVD, non-MUSCL; (b) Upwind TVD, non-MUSCL; (c) Upwind TVD, MUSCL; 
(d) van Leer splitting, MUSCL; (e) Steger and Warming splitting, MUSCL. 

Note: (a)-(c) use the identical approximate Riemann solver. 
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Comparison of the schemes for test case C, for equilibrium air. 

(a) Symmetric TVD, non-MUSCL; (b) Upwind TVD, non-MUSCL; (c) Upwind TVD, MUSCL, 
(d) Generalized van Leer splitting, MUSCL; (e) Generalized Steger and Warming splitting, MUSCL. 
Note: (a)-(c) use the identical approximate Riemann solver. 






Comparison of the schemes for test case D, for equilibrium air. 

(a) Symmetric TVD, non-MUSCL; (b) Upwind TVD, non-MUSCL; (c) Upwind TVD, MUSCL; 
(d) Generalized van Leer splitting, MUSCL; (e) Generalized Steger and Warming splitting, MUSCL 
Note: (a)-(c) use the identical approximate Riemann solver. 




Comparison of the schemes for test case D, for the perfect gas. 

(a) Symmetric TVD, non-MUSCL; (b) Upwind TVD, non-MUSCL; (c) Upwind TVD, MUSCL; 
(d) van Leer splitting, MUSCL; (e) Steger and Warming splitting, MUSCL. 

Note: (a)-(c) use the identical approximate Riemann solver. 
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Comparison of the schemes for test case E, for the perfect gas. 

(a) Symmetric TVD, non-MUSCL; (b) Upwind TVD, non-MUSCL; (c) Upwind TVD, MUSCL; 
(d) van Leer splitting, MUSCL; (e) Steger and Warming splitting, MUSCL. 

Note: (a)-(c) use the identical approximate Riemann solver. 







Fig. 13 Comparison of the schemes for test case F, for the perfect gas. 

(a) Symmetric TVD, non-MUSCL; (b) Upwind TVD, non-MUSCL; (c) Upwind TVD, MUSCL 
(d) van Leer splitting, MUSCL; (e) Steger and Warming splitting, MUSCL. 

Note: (a)-(c) use the identical approximate Riemann solver. 
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